﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;

namespace LowLevelGraphics
{
    /// <summary>
    /// Use this class for certain needed calculations
    /// </summary>
    public class ImageMath
    {
        /// <summary>
        /// The value of half pi as a float.
        /// </summary>
        public static float HALF_PI = (float)Math.PI / 2.0f;

        /// <summary>
        /// The value of quarter pi as a float.
        /// </summary>
        public static float QUARTER_PI = (float)Math.PI / 4.0f;

        /// <summary>
        /// The value of two pi as a float.
        /// </summary>
        public static float TWO_PI = (float)Math.PI * 2.0f;

        public struct Vector
        {
            public PointF Origin;
            public float Direction;

            /// <summary>
            /// Initializes a new instance of the <see cref="Vector"/> struct.
            /// </summary>
            /// <param name="_origin">The _origin.</param>
            /// <param name="_fDirection">The _f direction.</param>
            public Vector(PointF _origin, float _fDirection)
            {
                this.Origin = _origin;
                this.Direction = _fDirection;
            }
        }

        /// <summary>
        /// This is the equivalent to matlab diff
        /// </summary>
        /// <typeparam name="T"></typeparam>
        /// <param name="t"></param>
        /// <returns></returns>
        public int[] diff(int[] _aInt)
        {
            int nCount = _aInt.Length;
            int[] aDiff = new int[nCount - 1];

            for (int i = 0; i < nCount - 1; i++)
            {
                aDiff[i] = _aInt[i] - _aInt[i + 1];
            }
            return aDiff;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="_aFloat"></param>
        /// <returns></returns>
        public float[] diff(float[] _aFloat)
        {
            int nCount = _aFloat.Length;
            float[] aDiff = new float[nCount - 1];

            for (int i = 0; i < nCount - 1; i++)
            {
                aDiff[i] = _aFloat[i] - _aFloat[i + 1];
            }
            return aDiff;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="_aByte"></param>
        /// <returns></returns>
        public byte[] diff(byte[] _aByte)
        {
            int nCount = _aByte.Length;
            byte[] aDiff = new byte[nCount - 1];

            for (int i = 0; i < nCount - 1; i++)
            {
                aDiff[i] = (byte)(_aByte[i] - _aByte[i + 1]);
            }
            return aDiff;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="_aDouble"></param>
        /// <returns></returns>
        public double[] diff(double[] _aDouble)
        {
            int nCount = _aDouble.Length;
            double[] aDiff = new double[nCount - 1];

            for (int i = 0; i < nCount - 1; i++)
            {
                aDiff[i] = _aDouble[i] - _aDouble[i + 1];
            }
            return aDiff;
        }

        /// <summary>
        /// Cubics the interpolate.
        /// * Catmull-Rom cubic interpolation
        /// equally spaced points p0, p1, p2, p3
        /// <!-- interpolate 0<= u < 1 between p1 and p2 -->
        /// (1 u u^2 u^3)
        /// (  0.0  1.0  0.0  0.0 ) (p0)
        /// ( -0.5  0.0  0.5  0.0 ) (p1)
        /// (  1.0 -2.5  2.0 -0.5 ) (p2)
        /// ( -0.5  1.5 -1.5  0.5 ) (p3)
        /// </summary>
        /// <param name="src">The SRC.</param>
        /// <param name="nRowStride">The n row stride.</param>
        /// <param name="nSrcDepth">The n SRC depth.</param>
        /// <param name="_aDst">The _a DST.</param>
        /// <param name="_nDstDepth">The _n DST depth.</param>
        /// <param name="_dx">The _DX.</param>
        /// <param name="_dy">The _dy.</param>
        /// <param name="_dBrighten">The _d brighten.</param>
        public static void CubicInterpolate(byte[] src, int nRowStride, int nSrcDepth, byte[] _aDst,
                                           int _nDstDepth, double _dx, double _dy, double _dBrighten)
        {
            int MAX_PIXEL_DEPTH = 4;
            float um1 = 0.0f, u = 0.0f, up1 = 0.0f, up2 = 0.0f;
            float vm1 = 0.0f, v = 0.0f, vp1 = 0.0f, vp2 = 0.0f;
            int c = 0;
            float[] verts = new float[4 * MAX_PIXEL_DEPTH];

            um1 = (float)(((-0.5f * _dx + 1.0f) * _dx - 0.5f) * _dx);
            u = (float)((1.5f * _dx - 2.5f) * _dx * _dx + 1.0f);
            up1 = (float)(((-1.5f * _dx + 2.0f) * _dx + 0.5f) * _dx);
            up2 = (float)((0.5f * _dx - 0.5f) * _dx * _dx);

            vm1 = (float)(((-0.5f * _dy + 1.0f) * _dy - 0.5f) * _dy);
            v = (float)((1.5f * _dy - 2.5f) * _dy * _dy + 1.0f);
            vp1 = (float)(((-1.5f * _dy + 2.0f) * _dy + 0.5f) * _dy);
            vp2 = (float)((0.5f * _dy - 0.5f) * _dy * _dy);

            // Note: if dstDepth < srcDepth, we calculate unneeded pixels here 
            // later - select or create index array.

            for (c = 0; c < 4 * nSrcDepth; ++c)
            {
                verts[c] = vm1 * src[c] + v * src[c + nRowStride] + vp1 * src[c + nRowStride * 2] + vp2 * src[c + nRowStride * 3];
            }

            for (c = 0; c < _nDstDepth; ++c)
            {
                float result;
                result = um1 * verts[c] + u * verts[c + nSrcDepth] + up1 * verts[c + nSrcDepth * 2] + up2 * verts[c + nSrcDepth * 3];
                result = (float)(result * _dBrighten);

                if (result < 0.0)
                {
                    _aDst[c] = 0;
                }
                else if (result > 255.0)
                {
                    _aDst[c] = 255;
                }
                else
                {
                    _aDst[c] = (byte)result;
                }
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="x"></param>
        /// <param name="mean"></param>
        /// <param name="stddev"></param>
        /// <returns></returns>
        public static double GetDensityNormalDistribution(double x, double mean, double stddev)
        {
            double firstTerm = 1.0 / (stddev * Math.Sqrt(2.0 * Math.PI));
            double secondTerm = -((x - mean) * (x - mean) / (2.0 * stddev * stddev));
            double result = firstTerm * Math.Exp(secondTerm);
            return result;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public static double GetCumulativeDensityNormalDistribution(double x)
        {
            const double c0 = 0.2316419;
            const double c1 = 1.330274429;
            const double c2 = 1.821255978;
            const double c3 = 1.781477937;
            const double c4 = 0.356563782;
            const double c5 = 0.319381530;
            const double c6 = 0.398942280401;
            double negative = (x < 0 ? 1.0 : 0.0);
            double xPos = (x < 0.0 ? -x : x);
            double k = 1.0 / (1.0 + (c0 * xPos));
            double y1 = (((((((c1 * k - c2) * k) + c3) * k) - c4) * k) + c5) * k;
            double y2 = 1.0 - (c6 * Math.Exp(-0.5 * xPos * xPos) * y1);
            return ((1.0 - negative) * y2) + (negative * (1.0 - y2));
        }

        double GaborFunction(
              double dx,
              double dy,
              double angle,
              double frequency,
              double sigma)
        {
            //assert(sigma > 0.0);
            double gaborScaleX = Math.Cos(angle) * frequency; //Sometimes called 'k0x'
            double gaborScaleY = Math.Sin(angle) * frequency; //Sometimes called 'k0y'
            return (Math.Cos((gaborScaleX * dx) + (gaborScaleY * dy))

                                      * (1.0 / (Math.Sqrt(2.0 * Math.PI * sigma * sigma)))
                                      * Math.Exp(-1.0 * ((dx * dx) + (dy * dy)) / (2.0 * sigma * sigma))

              );
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="sigma"></param>
        /// <param name="min"></param>
        /// <param name="max"></param>
        public void GetMinMaxGaborFunction(double sigma, out double min, out double max)
        {
            //assert(sigma > 0.0);
            max = 1.0 / (Math.Sqrt(2.0 * Math.PI * sigma * sigma));
            min = -max;
        }

        /// <summary>
        /// Gets the distance.
        /// </summary>
        /// <param name="A">The A.</param>
        /// <param name="B">The B.</param>
        /// <returns></returns>
        public static double GetDistance(PointF A, PointF B)
        {
            float a = A.X - B.X;
            float b = A.Y - B.Y;
            return Math.Sqrt((double)(a * a + b * b));
        }

        private static double urand()
        {
            return urand(double.MaxValue);
        }

        private static double urand(double _dMaxRand)
        {
            return ((double)new Random().NextDouble() / _dMaxRand);
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="g1"></param>
        /// <param name="g2"></param>
        public static void GaussianRand(ref double g1, ref double g2)
        {
            GaussianRand(ref g1, ref g2, double.MaxValue);
        }

        /// <summary>
        /// Generate Gaussian Noise
        /// </summary>
        /// <param name="g1"></param>
        /// <param name="g2"></param>
        public static void GaussianRand(ref double g1, ref double g2, double _dMaxRand)
        {
            double M_PI = Math.PI;
            double n1 = urand(_dMaxRand);
            double n2 = urand(_dMaxRand);
            double x1 = n1 + ((n1 == 0) ? 1 : 0); /* guard against log(0) */
            double sqlogn1 = Math.Sqrt(-2.0 * Math.Log(x1));
            double angl = (2.0 * M_PI) * n2;
            g1 = sqlogn1 * Math.Cos(angl);
            g2 = sqlogn1 * Math.Sin(angl);
        }

        /// <summary>
        /// Classic euclidean algorithm.
        /// </summary>
        /// <param name="_x">The _x.</param>
        /// <param name="_y">The _y.</param>
        /// <returns></returns>
        public static int Euclidean(int _x, int _y)
        {
            while (_x != _y)
            {
                if (_x > _y)
                {
                    _x = _x - _y;
                }
                else
                {
                    _y = _y - _x;
                }
            }
            return _x;
        }

        /// <summary>
        /// Calculates euclidean.
        /// </summary>
        /// <param name="_x">The _x.</param>
        /// <param name="_y">The _y.</param>
        /// <returns></returns>
        public static int BinaryEuclidean(int _x, int _y)
        {
            int d = 1;
            while ((_x % 2 == 0) && (_y % 2 == 0))
            {
                _x = _x / 2; _y = _y / 2; d = 2 * d;
            }
            //and now the ggT algorithm
            while (_x != 0)
            {
                while (_x % 2 == 0) _x = _x / 2;
                while (_y % 2 == 0) _y = _y / 2;
                if (_x < _y) Swap<int>(ref _x, ref _y);
                _x = _x - _y;
            }
            return (_y * d);
        }

        /// <summary>
        /// Swaps the specified a.
        /// </summary>
        /// <typeparam name="T"></typeparam>
        /// <param name="_a">The _a.</param>
        /// <param name="_b">The _b.</param>
        public static void Swap<T>(ref T _a, ref T _b)
        {
            T temp = _a;
            _a = _b;
            _b = temp;
        }

        /// <summary>
        /// Gets the intersection.
        /// </summary>
        /// <param name="_vector1">The _vector1.</param>
        /// <param name="_vector2">The _vector2.</param>
        /// <returns></returns>
        public static PointF GetIntersection(Vector _vector1, Vector _vector2) //Vector[] pointAngularCoeff)
        {
            if (_vector1.Origin.X == _vector2.Origin.X && _vector1.Origin.Y == _vector2.Origin.Y)
                return _vector1.Origin;

            if (_vector1.Direction == _vector2.Direction) return PointF.Empty; //Parallel, no intersection

            float fNewX = float.Epsilon;
            float fNewY = float.Epsilon;

            if (float.IsInfinity(_vector1.Direction))
            {
                fNewX = _vector1.Origin.X;
                fNewY = _vector2.Origin.Y + _vector2.Direction * (-_vector2.Origin.X + _vector1.Origin.X);
            }

            if (float.IsInfinity(_vector2.Direction))
            {
                fNewX = _vector2.Origin.X;
                fNewY = _vector1.Origin.Y + _vector1.Direction * (-_vector1.Origin.X + _vector2.Origin.X);
            }

            if (fNewX == float.Epsilon)
            {
                float q1 = _vector1.Origin.Y - _vector1.Direction * _vector1.Origin.X;
                float q2 = _vector2.Origin.Y - _vector2.Direction * _vector2.Origin.X;
                fNewX = (q1 - q2) / (_vector2.Direction - _vector1.Direction);
                fNewY = _vector1.Direction * fNewX + q1;
            }

            if (float.IsInfinity(fNewX) || float.IsInfinity(fNewY))
                return PointF.Empty; //no intersection found
            else
            {
                return new PointF(fNewX, fNewY);
            }
        }

        /// <summary>
        /// Gets the angle.
        /// </summary>
        /// <param name="_from">The _from.</param>
        /// <param name="_to">The _to.</param>
        /// <returns></returns>
        public static float GetAngle(PointF _ptFrom, PointF _ptTo)
        {
            double angle = GetAngleRad(_ptFrom, _ptTo);
            double t = angle % Math.PI;

            if (t == HALF_PI)
                return float.PositiveInfinity;
            else
            {
                if (t == -HALF_PI)
                    return float.NegativeInfinity;
                else
                    return (float)Math.Tan(angle);
            }
        }

        /// <summary>
        /// Gets the angle RAD.
        /// </summary>
        /// <param name="_ptFrom">The _from.</param>
        /// <param name="_ptTo">The _to.</param>
        /// <returns></returns>
        public static double GetAngleRad(PointF _ptFrom, PointF _ptTo)
        {
            if (_ptTo.Y == _ptFrom.Y)
            {
                if (_ptFrom.X > _ptTo.X)
                    return Math.PI;
                else
                    return 0;
            }
            else if (_ptTo.X == _ptFrom.X)
            {
                if (_ptTo.Y < _ptFrom.Y)
                    return -HALF_PI;
                else
                    return HALF_PI;
            }
            else
            {
                double m = Math.Atan(((_ptTo.Y - _ptFrom.Y) / (_ptTo.X - _ptFrom.X)));

                if (_ptTo.X < 0)
                    if (m > 0)
                        return m + HALF_PI;
                    else
                        return m - Math.PI;
                else
                    return m;
            }
        }

        /// <summary>
        /// Calculates Linear interpolation.
        /// </summary>
        /// <param name="_t">The _t.</param>
        /// <param name="_a">The _a.</param>
        /// <param name="_b">The _b.</param>
        /// <returns></returns>
        public static float LinearInterpolate(float _t, float _a, float _b)
        {
            return _a + _t * (_b - _a);
        }

        /// <summary>
        /// Calculates Linear interpolation.
        /// </summary>
        /// <param name="_t">The _t.</param>
        /// <param name="_a">The _a.</param>
        /// <param name="_b">The _b.</param>
        /// <returns></returns>
        public static int LinearInterpolate(float _t, int _a, int _b)
        {
            return (int)(_a + _t * (_b - _a));
        }

        /// <summary>
        /// Clamps the specified c value.
        /// </summary>
        /// <param name="_c">The _c.</param>
        /// <returns></returns>
        public static int Clamp(int _c)
        {
            if (_c < 0)
                return 0;
            if (_c > 255)
                return 255;
            return _c;
        }

        /// <summary>
        /// Clamp a value to an interval.
        /// </summary>
        /// <param name="_x">The _x.</param>
        /// <param name="_a">The _a.</param>
        /// <param name="_b">The _b.</param>
        /// <returns>return the clamped value</returns>
        public static float Clamp(float _x, float _a, float _b)
        {
            return (_x < _a) ? _a : (_x > _b) ? _b : _x;
        }

        /// <summary>
        /// The triangle function. Returns a repeating triangle shape in the range 0..1 with wavelength 1.0
        /// </summary>
        /// <param name="x">the input parameter</param>
        /// <returns>the output value</returns>
        public static float Triangle(float x)
        {
            float r = Mod(x, 1.0f);
            return 2.0f * (r < 0.5 ? r : 1 - r);
        }

        /// <summary>
        /// Clamp a value to an interval.
        /// </summary>
        /// <param name="_x">The _x.</param>
        /// <param name="_a">The _a.</param>
        /// <param name="_b">The _b.</param>
        /// <returns>return the clamped value</returns>
        public static int Clamp(int _x, int _a, int _b)
        {
            return (_x < _a) ? _a : (_x > _b) ? _b : _x;
        }

        /// <summary>
        /// Apply a bias to a number in the unit interval, moving numbers towards 0 or 1
        /// according to the bias parameter.
        /// </summary>
        /// <param name="_fA">The _f A.</param>
        /// <param name="_fB">The _f B.</param>
        /// <returns>return the output value</returns>
        public static float Bias(float _fA, float _fB)
        {
            //		return (float)Math.pow(a, Math.log(b) / Math.log(0.5));
            return _fA / ((1.0f / _fB - 2) * (1.0f - _fA) + 1);
        }

        /// <summary>
        /// A variant of the gamma function..
        /// </summary>
        /// <param name="_fA">The _f A.</param>
        /// <param name="_fB">The _f B.</param>
        /// <returns>the output value</returns>
        public static float Gain(float _fA, float _fB)
        {
            /*
                    float p = (float)Math.log(1.0 - b) / (float)Math.log(0.5);

                    if (a < .001)
                        return 0.0f;
                    else if (a > .999)
                        return 1.0f;
                    if (a < 0.5)
                        return (float)Math.pow(2 * a, p) / 2;
                    else
                        return 1.0f - (float)Math.pow(2 * (1. - a), p) / 2;
            */
            float c = (1.0f / _fB - 2.0f) * (1.0f - 2.0f * _fA);
            if (_fA < 0.5)
                return _fA / (c + 1.0f);
            else
                return (c - _fA) / (c - 1.0f);
        }

        /// <summary>
        /// The pulse function. Returns 1 between two thresholds, 0 outside.
        /// </summary>
        /// <param name="_fA">The _f A.</param>
        /// <param name="_fB">The _f B.</param>
        /// <param name="_fX">The _f X.</param>
        /// <returns>the output value - 0 or 1</returns>
        public static float Pulse(float _fA, float _fB, float _fX)
        {
            return (_fX < _fA || _fX >= _fB) ? 0.0f : 1.0f;
        }

        /// <summary>
        /// A smoothed pulse function. A cubic function is used to smooth the step between two thresholds.
        /// </summary>
        /// <param name="a1">a1 the lower threshold position for the start of the pulse</param>
        /// <param name="a2">a2 the upper threshold position for the start of the pulse</param>
        /// <param name="b1">b1 the lower threshold position for the end of the pulse</param>
        /// <param name="b2">b2 the upper threshold position for the end of the pulse</param>
        /// <param name="x">x the input parameter</param>
        /// <returns>the output value</returns>
        public static float SmoothPulse(float a1, float a2, float b1, float b2, float x)
        {
            if (x < a1 || x >= b2)
                return 0;
            if (x >= a2)
            {
                if (x < b1)
                    return 1.0f;
                x = (x - b1) / (b2 - b1);
                return 1.0f - (x * x * (3.0f - 2.0f * x));
            }
            x = (x - a1) / (a2 - a1);
            return x * x * (3.0f - 2.0f * x);
        }

        /// <summary>
        /// The step function. Returns 0 below a threshold, 1 above.
        /// </summary>
        /// <param name="a">the threshold position</param>
        /// <param name="x">the input parameter</param>
        /// <returns>the output value - 0 or 1</returns>
        public static float Step(float a, float x)
        {
            return (x < a) ? 0.0f : 1.0f;
        }

        /// <summary>
        /// A smoothed step function. A cubic function is used to smooth the step between two thresholds.
        /// </summary>
        /// <param name="a">a the lower threshold position</param>
        /// <param name="b">b the upper threshold position</param>
        /// <param name="x">x the input parameter</param>
        /// <returns>the output value</returns>
        public static float SmoothStep(float a, float b, float x)
        {
            if (x < a)
                return 0;
            if (x >= b)
                return 1;
            x = (x - a) / (b - a);
            return x * x * (3 - 2 * x);
        }

        /// <summary>
        /// A "circle up" function. Returns y on a unit circle given 1-x. Useful for forming bevels.
        /// </summary>
        /// <param name="x">x the input parameter in the range 0..1</param>
        /// <returns> the output value</returns>
        public static float CircleUp(float x)
        {
            x = 1 - x;
            return (float)Math.Sqrt(1 - x * x);
        }

        /// <summary>
        /// A "circle down" function. Returns 1-y on a unit circle given x. Useful for forming bevels.
        /// </summary>
        /// <param name="x">the input parameter in the range 0..1</param>
        /// <returns>the output value</returns>
        public static float CircleDown(float x)
        {
            return 1.0f - (float)Math.Sqrt(1 - x * x);
        }

        /// <summary>
        /// Return a mod b. This differs from the % operator with respect to negative numbers.
        /// </summary>
        /// <param name="a">a the dividend</param>
        /// <param name="b">b the divisor</param>
        /// <returns> a mod b</returns>
        public static double Mod(double a, double b)
        {
            int n = (int)(a / b);

            a -= n * b;
            if (a < 0)
                return a + b;
            return a;
        }

        /// <summary>
        /// Return a mod b. This differs from the % operator with respect to negative numbers.
        /// </summary>
        /// <param name="a">a the dividend</param>
        /// <param name="b">b the divisor</param>
        /// <returns> a mod b</returns>
        public static float Mod(float a, float b)
        {
            int n = (int)(a / b);

            a -= n * b;
            if (a < 0)
                return a + b;
            return a;
        }

        /// <summary>
        /// Return a mod b. This differs from the % operator with respect to negative numbers.
        /// </summary>
        /// <param name="a">a the dividend</param>
        /// <param name="b">b the divisor</param>
        /// <returns> a mod b</returns>   
        public static int Mod(int a, int b)
        {
            int n = a / b;

            a -= n * b;
            if (a < 0)
                return a + b;
            return a;
        }


        // Catmull-Rom splines
        private static float m00 = -0.5f;
        private static float m01 = 1.5f;
        private static float m02 = -1.5f;
        private static float m03 = 0.5f;
        private static float m10 = 1.0f;
        private static float m11 = -2.5f;
        private static float m12 = 2.0f;
        private static float m13 = -0.5f;
        private static float m20 = -0.5f;
        private static float m21 = 0.0f;
        private static float m22 = 0.5f;
        private static float m23 = 0.0f;
        private static float m30 = 0.0f;
        private static float m31 = 1.0f;
        private static float m32 = 0.0f;
        private static float m33 = 0.0f;

        /// <summary>
        /// Compute a Catmull-Rom spline.
        /// </summary>
        /// <param name="x">x the input parameter</param>
        /// <param name="numKnots">numKnots the number of knots in the spline</param>
        /// <param name="_aKnots">The _a knots.</param>
        /// <returns>the spline value</returns>
        public static float Spline(float x, int numKnots, float[] _aKnots)
        {
            int span;
            int numSpans = numKnots - 3;
            float k0, k1, k2, k3;
            float c0, c1, c2, c3;

            if (numSpans < 1)
                throw new ArgumentException("Too few knots in spline");

            x = Clamp(x, 0, 1) * numSpans;
            span = (int)x;
            if (span > numKnots - 4)
                span = numKnots - 4;
            x -= span;

            k0 = _aKnots[span];
            k1 = _aKnots[span + 1];
            k2 = _aKnots[span + 2];
            k3 = _aKnots[span + 3];

            c3 = m00 * k0 + m01 * k1 + m02 * k2 + m03 * k3;
            c2 = m10 * k0 + m11 * k1 + m12 * k2 + m13 * k3;
            c1 = m20 * k0 + m21 * k1 + m22 * k2 + m23 * k3;
            c0 = m30 * k0 + m31 * k1 + m32 * k2 + m33 * k3;

            return ((c3 * x + c2) * x + c1) * x + c0;
        }

        /// <summary>
        /// Computes a Catmull-Rom spline, but with variable knot spacing.
        /// </summary>
        /// <param name="x">the input parameter</param>
        /// <param name="numKnots">the number of knots in the spline</param>
        /// <param name="_aXknots">The _a xknots.</param>
        /// <param name="_aYknots">The _a yknots.</param>
        /// <returns>the spline value</returns>
        public static float Spline(float x, int numKnots, int[] _aXknots, int[] _aYknots)
        {
            int span;
            int numSpans = numKnots - 3;
            float k0, k1, k2, k3;
            float c0, c1, c2, c3;

            if (numSpans < 1)
                throw new ArgumentException("Too few knots in spline");

            for (span = 0; span < numSpans; span++)
                if (_aXknots[span + 1] > x)
                    break;
            if (span > numKnots - 3)
                span = numKnots - 3;
            float t = (float)(x - _aXknots[span]) / (_aXknots[span + 1] - _aXknots[span]);
            span--;
            if (span < 0)
            {
                span = 0;
                t = 0;
            }

            k0 = _aYknots[span];
            k1 = _aYknots[span + 1];
            k2 = _aYknots[span + 2];
            k3 = _aYknots[span + 3];

            c3 = m00 * k0 + m01 * k1 + m02 * k2 + m03 * k3;
            c2 = m10 * k0 + m11 * k1 + m12 * k2 + m13 * k3;
            c1 = m20 * k0 + m21 * k1 + m22 * k2 + m23 * k3;
            c0 = m30 * k0 + m31 * k1 + m32 * k2 + m33 * k3;

            return ((c3 * t + c2) * t + c1) * t + c0;
        }

        /// <summary>
        /// Compute a Catmull-Rom spline for RGB values.
        /// </summary>
        /// <param name="x">the input parameter</param>
        /// <param name="numKnots">the number of knots in the spline</param>
        /// <param name="_aKnots">The _a knots.</param>
        /// <returns>the spline value</returns>
        public static int ColorSpline(float x, int numKnots, int[] _aKnots)
        {
            int span;
            int numSpans = numKnots - 3;
            float k0, k1, k2, k3;
            float c0, c1, c2, c3;

            if (numSpans < 1)
                throw new ArgumentException("Too few knots in spline");

            x = Clamp(x, 0, 1) * numSpans;
            span = (int)x;
            if (span > numKnots - 4)
                span = numKnots - 4;
            x -= span;

            int v = 0;
            for (int i = 0; i < 4; i++)
            {
                int shift = i * 8;

                k0 = (_aKnots[span] >> shift) & 0xff;
                k1 = (_aKnots[span + 1] >> shift) & 0xff;
                k2 = (_aKnots[span + 2] >> shift) & 0xff;
                k3 = (_aKnots[span + 3] >> shift) & 0xff;

                c3 = m00 * k0 + m01 * k1 + m02 * k2 + m03 * k3;
                c2 = m10 * k0 + m11 * k1 + m12 * k2 + m13 * k3;
                c1 = m20 * k0 + m21 * k1 + m22 * k2 + m23 * k3;
                c0 = m30 * k0 + m31 * k1 + m32 * k2 + m33 * k3;
                int n = (int)(((c3 * x + c2) * x + c1) * x + c0);
                if (n < 0)
                    n = 0;
                else if (n > 255)
                    n = 255;
                v |= n << shift;
            }

            return v;
        }

        /// <summary>
        /// Compute a Catmull-Rom spline for RGB values, but with variable knot space.
        /// </summary>
        /// <param name="x">the input parameter</param>
        /// <param name="numKnots">the number of knots in the spline</param>
        /// <param name="_aXknots">The _a xknots.</param>
        /// <param name="_aYknots">The _a yknots.</param>
        /// <returns>the spline value</returns>
        public static int ColorSpline(int x, int numKnots, int[] _aXknots, int[] _aYknots)
        {
            int span;
            int numSpans = numKnots - 3;
            float k0, k1, k2, k3;
            float c0, c1, c2, c3;

            if (numSpans < 1)
                throw new ArgumentException("Too few knots in spline");

            for (span = 0; span < numSpans; span++)
                if (_aXknots[span + 1] > x)
                    break;
            if (span > numKnots - 3)
                span = numKnots - 3;
            float t = (float)(x - _aXknots[span]) / (_aXknots[span + 1] - _aXknots[span]);
            span--;
            if (span < 0)
            {
                span = 0;
                t = 0;
            }

            int v = 0;
            for (int i = 0; i < 4; i++)
            {
                int shift = i * 8;

                k0 = (_aYknots[span] >> shift) & 0xff;
                k1 = (_aYknots[span + 1] >> shift) & 0xff;
                k2 = (_aYknots[span + 2] >> shift) & 0xff;
                k3 = (_aYknots[span + 3] >> shift) & 0xff;

                c3 = m00 * k0 + m01 * k1 + m02 * k2 + m03 * k3;
                c2 = m10 * k0 + m11 * k1 + m12 * k2 + m13 * k3;
                c1 = m20 * k0 + m21 * k1 + m22 * k2 + m23 * k3;
                c0 = m30 * k0 + m31 * k1 + m32 * k2 + m33 * k3;
                int n = (int)(((c3 * t + c2) * t + c1) * t + c0);
                if (n < 0)
                    n = 0;
                else if (n > 255)
                    n = 255;
                v |= n << shift;
            }

            return v;
        }


        /// <summary>
        /// Premultiplies the specified p. (Block of pixels)
        /// </summary>
        /// <param name="_aP">The _a P.</param>
        /// <param name="offset">The offset.</param>
        /// <param name="length">The length.</param>
        /// Premultiply a block of pixels
        public static void Premultiply(int[] _aP, int offset, int length)
        {
            length += offset;
            for (int i = offset; i < length; i++)
            {
                int rgb = _aP[i];
                int a = (rgb >> 24) & 0xff;
                int r = (rgb >> 16) & 0xff;
                int g = (rgb >> 8) & 0xff;
                int b = rgb & 0xff;
                float f = a * (1.0f / 255.0f);
                r = (int)(r * f);
                g = (int)(g * f);
                b = (int)(b * f);
                _aP[i] = (a << 24) | (r << 16) | (g << 8) | b;
            }
        }

        /// <summary>
        /// Unpremultiplies the specified p. (Block of pixels)
        /// </summary>
        /// <param name="_aP">The _a P.</param>
        /// <param name="offset">The offset.</param>
        /// <param name="length">The length.</param>
        public static void Unpremultiply(int[] _aP, int offset, int length)
        {
            length += offset;
            for (int i = offset; i < length; i++)
            {
                int rgb = _aP[i];
                int a = (rgb >> 24) & 0xff;
                int r = (rgb >> 16) & 0xff;
                int g = (rgb >> 8) & 0xff;
                int b = rgb & 0xff;
                if (a != 0 && a != 255)
                {
                    float f = 255.0f / a;
                    r = (int)(r * f);
                    g = (int)(g * f);
                    b = (int)(b * f);
                    if (r > 255)
                        r = 255;
                    if (g > 255)
                        g = 255;
                    if (b > 255)
                        b = 255;
                    _aP[i] = (a << 24) | (r << 16) | (g << 8) | b;
                }
            }
        }

        /// <summary>
        /// Resamples the specified source.
        /// An implementation of Fant's resampling algorithm.
        /// </summary>
        /// <param name="_aSource">The _a source.</param>
        /// <param name="_aDest">The _a dest.</param>
        /// <param name="length">the length of the scanline to resample</param>
        /// <param name="offset">the start offset into the arrays</param>
        /// <param name="stride">the offset between pixels in consecutive rows</param>
        /// <param name="out">an array of output positions for each pixel</param>
        public static void Resample(int[] _aSource, int[] _aDest, int length, int offset, int stride, float[] @out)
        {
            int i, j;
            float sizfac;
            float inSegment;
            float outSegment;
            int a, r, g, b, nextA, nextR, nextG, nextB;
            float aSum, rSum, gSum, bSum;
            float[] @in;
            int srcIndex = offset;
            int destIndex = offset;
            int lastIndex = _aSource.Length;
            int rgb;

            @in = new float[length + 2];
            i = 0;
            for (j = 0; j < length; j++)
            {
                while (@out[i + 1] < j)
                    i++;
                @in[j] = i + (float)(j - @out[i]) / (@out[i + 1] - @out[i]);
                //			in[j] = ImageMath.clamp( in[j], 0, length-1 );
            }
            @in[length] = length;
            @in[length + 1] = length;

            inSegment = 1.0f;
            outSegment = @in[1];
            sizfac = outSegment;
            aSum = rSum = gSum = bSum = 0.0f;
            rgb = _aSource[srcIndex];
            a = (rgb >> 24) & 0xff;
            r = (rgb >> 16) & 0xff;
            g = (rgb >> 8) & 0xff;
            b = rgb & 0xff;
            srcIndex += stride;
            rgb = _aSource[srcIndex];
            nextA = (rgb >> 24) & 0xff;
            nextR = (rgb >> 16) & 0xff;
            nextG = (rgb >> 8) & 0xff;
            nextB = rgb & 0xff;
            srcIndex += stride;
            i = 1;

            while (i <= length)
            {
                float aIntensity = inSegment * a + (1.0f - inSegment) * nextA;
                float rIntensity = inSegment * r + (1.0f - inSegment) * nextR;
                float gIntensity = inSegment * g + (1.0f - inSegment) * nextG;
                float bIntensity = inSegment * b + (1.0f - inSegment) * nextB;
                if (inSegment < outSegment)
                {
                    aSum += (aIntensity * inSegment);
                    rSum += (rIntensity * inSegment);
                    gSum += (gIntensity * inSegment);
                    bSum += (bIntensity * inSegment);
                    outSegment -= inSegment;
                    inSegment = 1.0f;
                    a = nextA;
                    r = nextR;
                    g = nextG;
                    b = nextB;
                    if (srcIndex < lastIndex)
                        rgb = _aSource[srcIndex];
                    nextA = (rgb >> 24) & 0xff;
                    nextR = (rgb >> 16) & 0xff;
                    nextG = (rgb >> 8) & 0xff;
                    nextB = rgb & 0xff;
                    srcIndex += stride;
                }
                else
                {
                    aSum += (aIntensity * outSegment);
                    rSum += (rIntensity * outSegment);
                    gSum += (gIntensity * outSegment);
                    bSum += (bIntensity * outSegment);
                    _aDest[destIndex] =
                        ((int)Math.Min(aSum / sizfac, 255) << 24) |
                        ((int)Math.Min(rSum / sizfac, 255) << 16) |
                        ((int)Math.Min(gSum / sizfac, 255) << 8) |
                        (int)Math.Min(bSum / sizfac, 255);
                    destIndex += stride;
                    aSum = rSum = gSum = bSum = 0.0f;
                    inSegment -= outSegment;
                    outSegment = @in[i + 1] - @in[i];
                    sizfac = outSegment;
                    i++;
                }
            }
        }
    }
}
